ndarrayのfortran flagと実際のメモリ配列について

ndarrayのfortran flagと実際のメモリ配列について

要約

  • ndarrayで.Tで転置を行った場合、実際のメモリ上の配列は変更されず、.flagsで確認できるC形式配列(row-major)であるかFortran形式配列(column-major)であるかのフラグが入れ替わるのみ
  • np.saveで保存する場合も上記フラグは保たれたまま保存されるため、実際のレイアウトを転置させたい場合はndarrayに対し.copy()を読んだ上で保存する必要がある

実験

In [1]:
import numpy as np

2x2のndarrayを生成する。

In [2]:
a = np.array([[1,2], [3,4]])
a
Out[2]:
array([[1, 2],
       [3, 4]])

.Tで転置を取ると行と列が入れ替わる。

In [3]:
a.T
Out[3]:
array([[1, 3],
       [2, 4]])

元のndarrayの.flagsを表示すると、C_CONTIGUOUSTrueとなっている一方、F_CONTIGUOUSがFalseとなっていることが分かる。C_CONTIGUOUSはC形式配列であるかを表すフラグであり、F_CONTIGUOUSはFortran形式配列であるかを表すフラグである。

In [4]:
a.flags
Out[4]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

一方.Tを呼んだ後のndarrayの.flagsを表示すると、C_CONTIGUOUSFalseに変化し、一方F_CONTIGUOUSがTrueとなっている。

In [5]:
a.T.flags
Out[5]:
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

np.saveでの保存時の形式

.Tで転置を取ると.flagsのC形式配列であるかFortran形式配列であるかのフラグが入れ替わることが分かったが、実際の格納形式がどうなっているかを知るために、np.saveを用いて保存した.npyファイルに対して!hexdump -Cvを実行することでバイナリの中身を調べる。ここでhexdumpコマンドの-CはASCII表示を行うためのオプションであり、-vは全データを表示するためのオプションである(おそらくBSD系のhexdumpのみ)。

まず、np.saveで転置前後のndarrayを保存する。

In [6]:
np.save("a.npy", a)
np.save("a.T.npy", a.T)

これらをnp.loadで再度読み込み直すと、.flagsは読み込んだ後も保持されていることが分かる。

In [7]:
np.load("a.npy").flags
Out[7]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False
In [8]:
np.load("a.T.npy").flags
Out[8]:
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

ディスクへダンプした際の実際の格納形式がどうなっているかをhexdumpコマンドを用いて調べる。

In [9]:
!hexdump -Cv a.npy
00000000  93 4e 55 4d 50 59 01 00  76 00 7b 27 64 65 73 63  |.NUMPY..v.{'desc|
00000010  72 27 3a 20 27 3c 69 38  27 2c 20 27 66 6f 72 74  |r': '<i8', 'fort|
00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 46 61 6c 73  |ran_order': Fals|
00000030  65 2c 20 27 73 68 61 70  65 27 3a 20 28 32 2c 20  |e, 'shape': (2, |
00000040  32 29 2c 20 7d 20 20 20  20 20 20 20 20 20 20 20  |2), }           |
00000050  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
00000060  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
00000070  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 0a  |               .|
00000080  01 00 00 00 00 00 00 00  02 00 00 00 00 00 00 00  |................|
00000090  03 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
000000a0
In [10]:
!hexdump -Cv a.T.npy
00000000  93 4e 55 4d 50 59 01 00  76 00 7b 27 64 65 73 63  |.NUMPY..v.{'desc|
00000010  72 27 3a 20 27 3c 69 38  27 2c 20 27 66 6f 72 74  |r': '<i8', 'fort|
00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 54 72 75 65  |ran_order': True|
00000030  2c 20 27 73 68 61 70 65  27 3a 20 28 32 2c 20 32  |, 'shape': (2, 2|
00000040  29 2c 20 7d 20 20 20 20  20 20 20 20 20 20 20 20  |), }            |
00000050  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
00000060  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
00000070  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 0a  |               .|
00000080  01 00 00 00 00 00 00 00  02 00 00 00 00 00 00 00  |................|
00000090  03 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
000000a0

ヘッダーのfortran_orderが変化していることが分かる。更にbashのプロセス置換で両者のdiffを見てみる。

In [11]:
%%bash
diff <(hexdump -Cv a.npy) <(hexdump -Cv a.T.npy)
3,5c3,5
< 00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 46 61 6c 73  |ran_order': Fals|
< 00000030  65 2c 20 27 73 68 61 70  65 27 3a 20 28 32 2c 20  |e, 'shape': (2, |
< 00000040  32 29 2c 20 7d 20 20 20  20 20 20 20 20 20 20 20  |2), }           |
---
> 00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 54 72 75 65  |ran_order': True|
> 00000030  2c 20 27 73 68 61 70 65  27 3a 20 28 32 2c 20 32  |, 'shape': (2, 2|
> 00000040  29 2c 20 7d 20 20 20 20  20 20 20 20 20 20 20 20  |), }            |

ヘッダーのみ差分が存在し、実際のデータが格納されている後半部は変化していないことが分かる。

転置した状態のまま保存する方法

ndarrayに対して.copy()を呼んだ上で保存すれば良い。下記に実例を示す。 まず、.copy()を呼んだ後の.flagsを表示する。

In [12]:
a.T.copy().flags
Out[12]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

転置したにもかかわらずC_CONTIGUOUSがTrueになっていることが分かる。 更にこのndarrayをnp.saveで保存し、hexdumpで中身を見てみる。

In [13]:
np.save("a.T.copy.npy", a.T.copy())
In [14]:
!hexdump -C a.T.copy.npy
00000000  93 4e 55 4d 50 59 01 00  76 00 7b 27 64 65 73 63  |.NUMPY..v.{'desc|
00000010  72 27 3a 20 27 3c 69 38  27 2c 20 27 66 6f 72 74  |r': '<i8', 'fort|
00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 46 61 6c 73  |ran_order': Fals|
00000030  65 2c 20 27 73 68 61 70  65 27 3a 20 28 32 2c 20  |e, 'shape': (2, |
00000040  32 29 2c 20 7d 20 20 20  20 20 20 20 20 20 20 20  |2), }           |
00000050  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 20  |                |
*
00000070  20 20 20 20 20 20 20 20  20 20 20 20 20 20 20 0a  |               .|
00000080  01 00 00 00 00 00 00 00  03 00 00 00 00 00 00 00  |................|
00000090  02 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
000000a0

C_CONTIGUOUSがTrue、F_CONTIGUOUSがFalseになっていたことからも分かるように、ヘッダーのfortran_orderがFalseになっている。

当然、再度np.loadで読み込んだ場合でもこの.flagsは保たれている。

In [15]:
np.load("a.T.copy.npy").flags
Out[15]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

元のndarrayと、転置した上で.copy()を呼んだ後のndarrayをnp.saveしたバイナリを比較してみる。

In [16]:
%%bash
diff <(hexdump -Cv a.npy) <(hexdump -Cv a.T.copy.npy)
9,10c9,10
< 00000080  01 00 00 00 00 00 00 00  02 00 00 00 00 00 00 00  |................|
< 00000090  03 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
---
> 00000080  01 00 00 00 00 00 00 00  03 00 00 00 00 00 00 00  |................|
> 00000090  02 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|

両者ともにC形式配列となっているため、純粋にデータ部の並びのみが異なる。

同様に転置を取った場合と、その上で.copy()を呼んだ場合のバイナリを比較してみる。

In [17]:
%%bash
diff <(hexdump -Cv a.T.npy) <(hexdump -Cv a.T.copy.npy)
3,5c3,5
< 00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 54 72 75 65  |ran_order': True|
< 00000030  2c 20 27 73 68 61 70 65  27 3a 20 28 32 2c 20 32  |, 'shape': (2, 2|
< 00000040  29 2c 20 7d 20 20 20 20  20 20 20 20 20 20 20 20  |), }            |
---
> 00000020  72 61 6e 5f 6f 72 64 65  72 27 3a 20 46 61 6c 73  |ran_order': Fals|
> 00000030  65 2c 20 27 73 68 61 70  65 27 3a 20 28 32 2c 20  |e, 'shape': (2, |
> 00000040  32 29 2c 20 7d 20 20 20  20 20 20 20 20 20 20 20  |2), }           |
9,10c9,10
< 00000080  01 00 00 00 00 00 00 00  02 00 00 00 00 00 00 00  |................|
< 00000090  03 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|
---
> 00000080  01 00 00 00 00 00 00 00  03 00 00 00 00 00 00 00  |................|
> 00000090  02 00 00 00 00 00 00 00  04 00 00 00 00 00 00 00  |................|

この場合、ヘッダーもデータ部も異なっていることが分かる。

なぜこの記事を書いたか

とある.npyファイルのデータ部分を直接読み込むCプログラムを触っていた際に、.Tで転置した.npyを読み込んでも読み込んだデータ部分が全く変化していないということが発生したためである。 当該プログラムではヘッダーのfortran_orderを完全に無視していたため、.Tしただけのndarrayを保存した.npyではデータ部分が転置されておらず、意図しない動作となってしまっていた。これを解決するためには前述のように.copy()を呼んだ上で保存した.npyファイルを読み込む必要があった。

Travis CIによるNikolaブログ構築の自動化

Travis CIによるNikolaブログ構築の自動化

Nikolaでブログを構築するための方法は過去の記事(1, 2, 3)に書いていますが、手元でのビルドのためにpythonやnikolaがインストールされたPCが必要になるため、ブログ記事を書くための環境が限定されてしまうという問題がありました。 そこで、この記事ではNikola公式の記事を参考に、Travis CIを用いることでsrcブランチにブログ記事をコミットするだけでTravis CI側で自動的にブログを構築し、masterブランチにプッシュしてくれる仕組みを構築する過程を紹介します。 なお、以下の作業はすべてsrcブランチ上で行います。

conf.pyの編集

まず、nikola github_deployコマンド実行時に、デプロイと同時にsrcブランチもコミットするかどうかを制御するためのオプションをFalseに設定します。

GITHUB_COMMIT_SOURCE = False

これはTravis CIではソースブランチへのコミットをトリガーにnikola build && nikola github_deployを実行するため、nikola github_deployによってsrcブランチがコミットされると再度それがトリガーとなって無限ループに陥ってしまうことを防ぐためです。

.travis.ymlの作成

nikolaブログのルートディレクトリ(conf.pyが置いてあるディレクトリ)に以下のような.travis.ymlファイルを作成します。 これはNikola公式の記事に記載のtravis.ymlを改変したものになります。

language: python
cache: apt
sudo: false
addons:
  apt:
    packages:
    - language-pack-ja-base
    - language-pack-ja
branches:
  only:
  - src
python:
- 3.6
before_install:
- git config --global user.name 'USERNAME'
- git config --global user.email 'travis@invalid'
- git config --global push.default 'simple'
- pip install --upgrade pip wheel
- echo -e 'Host github.com\n    StrictHostKeyChecking no' >> ~/.ssh/config
- eval "$(ssh-agent -s)"
- chmod 600 id_rsa
- ssh-add id_rsa
- git remote rm origin
- git remote add origin git@github.com:USERNAME/REPO.git
- git fetch origin master
- git branch master FETCH_HEAD
install:
- pip install 'ghp-import2'
- pip install 'webassets'
- pip install 'Nikola[extras]'
script:
- nikola build && nikola github_deploy -m 'Nikola auto deploy [ci skip]'
notifications:
  email:
    on_success: change
    on_failure: always

ここで、

- git config --global user.name 'USERNAME'
- git config --global user.email 'travis@invalid'

の行は適切なユーザー名、メールアドレスに、

- git remote add origin git@github.com:USERNAME/REPO.git

の行は適切なユーザー名およびリポジトリ名に変更する必要があります。

SSH鍵の生成

まず、.gitignoreファイルにid_rsaid_rsa.pubを無視する設定を追記した上でSSH鍵を生成します。

$ echo id_rsa >> .gitignore
$ echo id_rsa.pub >> .gitignore
$ ssh-keygen -C TravisCI -f id_rsa -N ''

上記を実行すると秘密鍵id_rsaおよび公開鍵id_rsa.pubが作成されます。 念のため、.gitignoreの無視設定が合っているかを確認するためにgit statusを実行してもaddされる対象となっていないことを確認しておきましょう。 なおssh-keygen-Cはコメント、-fは鍵名、-Nはパスフレーズの指定(ここでは空文字なので指定なし)です。

公開鍵のgithubリポジトリへの登録

生成した公開鍵はgithubリポジトリに登録しておく必要があります。 リポジトリページ -> Settings -> Deploy Keys -> Add deploy key からTitleをTravis CIとし、Keyにidrsa.pubの中身をコピペして登録しましょう。 また、Allow write accessはチェックしておく必要があります。 これらの作業を忘れるとTravis CIの自動ビルド時にアクセスエラーが発生します。

travis encrypt-fileコマンドによる秘密鍵の暗号化

生成した秘密鍵id_rsatravis encrypt-fileによってid_rsa.encに共通鍵による暗号化を施した上でリポジトリに追加します。 この暗号化されたid_rsa.encはTravis CIでの自動ビルド時に共通鍵によって復号化され、Travis CIのgithubへのアクセスに使用されます。

この作業を実行するためにはtravisコマンドが必要となりますが、これはrubyのgemとして配布されているため、インストールされていない場合は下記にようにgemでインストールする必要があります。

$ gem install --user-install travis

なおtravisコマンドがインストールされる場所にPATHが通っていない場合はPATHに追加するか、下記一連のコマンドをフルパスで実行する必要があります。

travisコマンドが実行可能になったら、更に下記を実行します。

$ travis login
$ travis enable
$ travis encrypt-file id_rsa --add

travis loginを実行するとgithubアカウントのユーザ名、パスワードを求められるため入力してください。 この上でtravis enableを実行すると自動ビルドを有効化するリポジトリが正しいか確認されるため、正しければyesと入力しましょう。 更にtravis encrypt-file id_rsa --addを実行すると、秘密鍵id_rsaが暗号化されてid_rsa.encが生成されます。 このid_rsa.enc.gitignoreに追加されていないため、git addによってgitの管理下に置かれることになります。 なお引数の--addを付けておくと、.travis.ymlファイルにid_rsa.encの復号化を行うための下記のようなopensslのコマンドを追加してくれます。

before_install:
- openssl aes-256-cbc -K $encrypted_XXXXXX_key -iv $encrypted_XXXXXX_iv
  -in id_rsa.enc -out id_rsa -d

-Kは共通鍵、-ivは初期ベクトルの指定であり、指定されている値はTravis CI側で環境変数として設定されています。 (ブラウザからTravis CIの設定を見ると確認することができます)

srcブランチへの各種ファイルのadd, commitおよびgithubへのpush

以上の作業によりconf.py.gitignore.travis.ymlid_rsa.encの4つのファイルが編集・生成されているため、これをsrcブランチにaddした上でcommitし、更にremoteとなっているgithubにpushします。

$ git add .
$ git commit -am "Automate builds with Travis CI"
$ git push origin src

これによりTravis CIの自動ビルドが実行されるはずなので、あとはブラウザからTravis CIのページを確認し、ビルドが通っているかを確認するのが良いでしょう。

上記設定以降の記事の追加方法

上記までで設定した方法によってsrcブランチが変更される度にTravis CIがnikola buildおよびnikola github_deployを行ってくれるようになったため、postsディレクトリ以下に新規記事を追加したら、あとはこの新規記事をsrcブランチにgit add posts/XX.mdのようにgit addした上でgit commitして、更にgithubにgit push origin srcすれば自動的にブログがビルドされます。

参考

2018/5/19にPycon Osakaで'librosaで始める音楽情報処理'というタイトルで発表しました

2018/5/19にPyCon mini Osakaで「librosaで始める音楽情報検索」というタイトルで発表しました

公開が遅くなりましたが先月開催されたPyCon mini Osakaで発表した内容を公開します。この記事もAzure Notebooksで既に公開済ではありましたが、nikolaのipynb公開機能を使ってブログポストとして投稿しておきたいと思います。

librosaで始める音楽情報検索

  • 大橋 宏正(twitter: @wrist)

この資料

自己紹介

  • 大橋宏正(@wrist)
    • 某メーカで音響信号処理の業務に従事
    • 2013年ごろから大阪PRML読書会を主催(現在休止中)
    • PyData Osakaオーガナイザ

PyCon mini Osakaのテーマ

  • 「Pythonでなにしよう?」
    • Webアプリ
    • 科学技術計算
  • Pythonで音楽の分析を行いたい!

この発表

  • Pythonで音楽情報検索(MIR)を行うためのライブラリlibrosaの紹介

音楽情報検索(MIR)とは?

  • Music Information Retrievalの略
  • 音楽情報に関する処理全般
    • ビートトラッキング、楽器音分離、音楽推薦など
  • コミュニティ
    • ismir.net
    • 音楽情報科学研究会(www.sigmus.jp)

MIRをPythonで行う

  • Pythonで音楽データの処理が必要
  • 読み書き
    • wave, audioread, pysoundfile, scipy.io, ...
  • 線形代数、信号処理、フーリエ変換
    • numpy, scipy(signal, fftpack, ...)
  • 可視化
    • matplotlib, seaborn
  • これらを組み合わせることでさまざまな処理を実現

参考

音響信号処理についてはQiitaに過去に諸々記事を書いてます

librosa

librosaの特徴

  • MIRに必要となる典型的な処理が多く実装
  • 低い学習障壁
  • 一貫性、後方互換性、モジュール性
In [1]:
# pipでインストール
!pip install librosa
Requirement already satisfied: librosa in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (0.6.0)
Requirement already satisfied: scipy>=0.14.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.1.0)
Requirement already satisfied: resampy>=0.2.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.2.0)
Requirement already satisfied: numpy>=1.8.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.14.2)
Requirement already satisfied: decorator>=3.0.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (4.3.0)
Requirement already satisfied: audioread>=2.0.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (2.1.5)
Requirement already satisfied: six>=1.3 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.11.0)
Requirement already satisfied: scikit-learn!=0.19.0,>=0.14.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.19.1)
Requirement already satisfied: joblib>=0.7.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.11)
Requirement already satisfied: numba>=0.32 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from resampy>=0.2.0->librosa) (0.38.0)
Requirement already satisfied: llvmlite>=0.23.0dev0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from numba>=0.32->resampy>=0.2.0->librosa) (0.23.0)

librosaの具体的な紹介

  • 適宜デモを交えながらlibrosaを紹介
  • 発表の流れ
    • 音データの短時間周波数分析
    • 特徴量抽出
    • オンセット検出、ビートトラッキング
In [2]:
# 本発表では以下を仮定
import numpy as np
import scipy as sp
import scipy.signal as sg
import scipy.fftpack as fft
import matplotlib as mpl
import matplotlib.pyplot as plt
import IPython.display
In [3]:
# version
print("numpy version", np.version.full_version)
print("scipy version", sp.version.full_version)
print("matplotlib version", mpl.__version__)
numpy version 1.14.2
scipy version 1.1.0
matplotlib version 2.2.2
In [4]:
# librosa以下の各モジュールをimport
import librosa
import librosa.display
import librosa.feature
import librosa.effects
import librosa.onset
import librosa.beat

サウンドデータの読み込み

  • librosa.load
    • ファイルパスを与えることでデータを読み込み
    • デフォルトでモノラル化、サンプリングレート22050Hzとなる
    • バックエンドにaudioreadライブラリを使用
  • librosa.display.waveplotで波形データの可視化
In [5]:
# デモ音源のロード&再生
# y: モノラル音声データ, sr: サンプリングレート
y, sr = librosa.load(librosa.util.example_audio_file())
IPython.display.Audio(data=y, rate=sr)
Out[5]:
In [6]:
# 波形の描画
# 縦軸レンジを自動的に指定し、横軸を秒数表示にしてくれる
librosa.display.waveplot(y, sr)
Out[6]:
<matplotlib.collections.PolyCollection at 0x115239668>

サウンドデータのデジタル表現

  • 現実の音はアナログデータ
    • 時間方向にも振幅方向にも連続
  • コンピュータで音を取り扱うためには離散化が必要
    • 時間方向: サンプリング
    • 振幅方向: 量子化
In [7]:
b_pt = int(10.0*sr); e_pt = int((10.0+0.002)*sr)
fig = plt.figure(1, figsize=(14,6)); ax = fig.add_subplot(111)
ax.plot(y[b_pt:e_pt]); ax.plot(y[b_pt:e_pt], "o"); ax.grid();
ax.set_xlabel("time [pt]"); ax.set_ylabel("amplitude")
Out[7]:
Text(0,0.5,'amplitude')
In [8]:
# サンプリングレート: 1秒辺りのサンプル数
print("sampling rate: {0}".format(sr))
print("length: {0}[pt]={1}[s]".format(y.shape[0], y.shape[0]/sr))
sampling rate: 22050
length: 1355168[pt]=61.45886621315193[s]

音声データの周波数分析

  • 音声データの長さはデータにより異なる
  • 固定長の短時間波形に切り出した上で処理を行う
  • 短時間波形に対し離散フーリエ変換を適用することで周波数分析を行う
In [9]:
# strideトリックを使った信号を短時間波形フレームの配列に変換する関数
from numpy.lib.stride_tricks import as_strided

def split_frames(xs, frame_width=1024, noverlap=2):

    nsample = np.alen(xs)
    nshift = (frame_width // noverlap)
    nframe = int(nsample / nshift) - noverlap + 1
    
    nbyte = xs.dtype.itemsize

    # as_strided(xs, ys_shape, new_stride)
    xs_split = as_strided(xs, (nframe, frame_width), (nbyte * nshift, nbyte))
    
    return xs_split
In [10]:
def plot_overlap_waves():
    import matplotlib.gridspec as gridspec

    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    gs = gridspec.GridSpec(3, nframes + 1)
    
    ax_wav = fig.add_subplot(gs[0, :])
    ax_wav.plot(np.arange(half_width*(nframes+1))/sr, 
                y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)/sr])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    win = sg.get_window("hanning", int(frame_width))
    xs = np.arange(frame_width)

    for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:nframes]):
        y_win = win * y_cut
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width)/sr, 0.2 * win)
        ax = None
        if i % 2 == 0:
            ax = fig.add_subplot(gs[1, i:i+2])
            ax.plot(xs[:half_width], y_win[:half_width])
            ax.plot(xs[half_width:], y_win[half_width:])
            ax.set_xlim([0, frame_width])
        else:
            ax = fig.add_subplot(gs[2, i:i+2])
            ax.plot(xs[half_width:], y_win[half_width:])
            ax.plot(xs[:half_width], y_win[:half_width])
            ax.set_xlim([0, frame_width])

        ax.tick_params(labelleft="off",left="off") # y軸の削除
        ax.set_ylim([-ymax, ymax]); ax.grid()

    fig.subplots_adjust(wspace=0.0)
In [11]:
# 周波数分析時には通常窓関数を乗算した短時間波形を使用
# 切り出しはここでは1/2区間をオーバーラップさせている
plot_overlap_waves()
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

用語の整理

  • sr = sampling rate [Hz = pt/s]
    • 22050 [Hz]
  • frame = 各短時間波形
  • n_fft = 各フレームに含まれるサンプル数
    • 2048 [pt]
  • hop_length = フレーム間のサンプル数
    • 512 [pt]
In [12]:
def plot_overlap_freqs():
    import matplotlib.gridspec as gridspec

    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    gs = gridspec.GridSpec(3, nframes + 1)
    
    ax_wav = fig.add_subplot(gs[0, :])
    ax_wav.plot(y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    win = sg.get_window("hanning", int(frame_width))
    xs = np.arange(frame_width)
    w = np.linspace(0.0, 2.0 * np.pi, frame_width, endpoint=False)
    freqs = w[:frame_width//2+1] * (sr / 2) / np.pi
    for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
        y_win = win * y_cut
        Y_abs = np.abs(fft.fft(y_win, frame_width)[:frame_width//2+1])
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width), 0.2 * win)
        ax = None
        if i % 2 == 0:
            ax = fig.add_subplot(gs[1, i:i+2])
            ax.semilogy(Y_abs, freqs)
        else:
            ax = fig.add_subplot(gs[2, i:i+2])
            ax.semilogy(Y_abs, freqs)

        if i > 1:
            ax.tick_params(labelleft="off",left="off") # y軸の削除
        else:
            ax.tick_params(labelbottom="off",bottom="off") # x軸の削除
            ax.set_ylabel("frequency [Hz]")

        ax.set_xlabel("amplitude")
        ax.set_ylim([10, sr/2])
        ax.grid()

    fig.subplots_adjust(wspace=0.0)
In [13]:
# 各短時間波形を周波数分析したグラフを描画
# 縦軸が周波数、横軸が振幅スペクトル
plot_overlap_freqs()
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [14]:
def plot_spectrogram():
    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    
    ax_wav = fig.add_subplot(211)
    ax_wav.plot(np.arange(half_width*(nframes+1))/sr,
                y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)/sr])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    win = sg.get_window("hanning", int(frame_width))
    xs = np.arange(frame_width)
    w = np.linspace(0.0, 2.0 * np.pi, frame_width, endpoint=False)
    freqs = w[:frame_width//2+1] * (sr / 2) / np.pi

    Zs = []
    for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
        y_win = win * y_cut
        Y_abs = np.abs(fft.fft(y_win, frame_width)[:frame_width//2+1])
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width), 0.2 * win)
        Zs.append(20.0*np.log10(Y_abs))
    
    ax = fig.add_subplot(212)
    Xs = np.arange((nframes+1)) / sr / 2
    pcol = ax.pcolormesh(Xs, freqs, np.array(Zs).T, cmap='magma')
    cax = fig.colorbar(pcol)
    cax.set_clim(-80, 20)

    fig.subplots_adjust(wspace=0.0)
In [15]:
# 短時間波形を周波数分析したものをスペクトログラムとして表示
plot_spectrogram()

librosaにおける短時間周波数分析

  • 短時間離散フーリエ変換(Short time fourie transform)
    • librosa.stft
  • Constant-Q Transform(CQT)
    • librosa.cqt
  • 分析結果の可視化
    • librosa.specshow
In [16]:
# librosaでは短時間波形分析をメソッド一発で実行可能
D = librosa.stft(y)
In [17]:
# D[k, l]
# k: 周波数ビン  l: フレームインデックス
# kは実信号に対する2048ptフーリエ変換の結果なので1025ptのみ保持
D.shape
Out[17]:
(1025, 2647)

離散フーリエ変換$D[k,l]$の算出法

$$D[k, l]=\sum_{n=0}^{{\textrm N_{\rm fft}-1}} w[n]\cdot x\left[\left(l\cdot {\rm M}-\frac{N_{\rm fft}}{2}\right)+n\right]\exp\left\{-j\frac{2\pi k n}{N_{\rm fft}}\right\}$$

* $x[n]$は入力信号、$w[n]$は窓関数
* 周波数ビン$k$、フレームインデックス$l$、$M$はhop length
* $lM$を中心に前後$N/2$点に渡り離散フーリエ変換を実行
In [18]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_power = librosa.amplitude_to_db(np.abs(D), ref=np.max)
librosa.display.specshow(log_power, x_axis="time", y_axis="linear")
plt.colorbar()
Out[18]:
<matplotlib.colorbar.Colorbar at 0x115a2aba8>
In [19]:
# 縦軸を対数表示
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_power = librosa.amplitude_to_db(np.abs(D), ref=np.max)
librosa.display.specshow(log_power, x_axis="time", y_axis="log")
plt.colorbar()
Out[19]:
<matplotlib.colorbar.Colorbar at 0x1c3344cef0>

Constant-Q変換

$$C[k]= \frac{1}{N[k]} \sum_{n=0}^{{\textrm N[k]-1}} w[n, k]\cdot x\left[n\right]\exp\left\{-j\frac{2\pi Q[k] n}{N[k]}\right\}$$

  • 特徴
    • 分析する周波数間隔を対数的に配置
    • 低域は細かく分析し、高域は荒く分析
    • 音楽信号の分析に適する
In [20]:
C = librosa.cqt(y, sr)
In [21]:
# STFTと比較して縦軸の周波数データ数は圧倒的に少ない
# STFTは周波数間隔が等幅なのに対してCQTだと対数間隔となるため
print(D.shape)
print(C.shape)
(1025, 2647)
(84, 2647)
In [22]:
fmin = 32.70  # C1に相当するHz
fmax = fmin * ((2**(1/12)) ** 84)
print(fmax)
print(librosa.hz_to_note(fmin))
print(librosa.hz_to_note(fmax))
4185.600000000015
C1
C8
In [23]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
plt.colorbar()
Out[23]:
<matplotlib.colorbar.Colorbar at 0x115933cf8>
In [24]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_note")
plt.colorbar()
Out[24]:
<matplotlib.colorbar.Colorbar at 0x11589d518>
In [25]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

# デフォルトだと80dBまでだがtop_dbを付けるとレンジが縮まる
log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max, top_db=40)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_note")
plt.colorbar()
Out[25]:
<matplotlib.colorbar.Colorbar at 0x1c1f823518>

クロマベクトル

  • 12音階それぞれの強さを要素に持つベクトル
  • CQTで算出した対数周波数帯域ごとの強度を元に計算
  • librosa.feature.chroma_cqt
In [26]:
chroma = librosa.feature.chroma_cqt(C=C, sr=sr)

fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

librosa.display.specshow(chroma, x_axis="time", y_axis="chroma")
plt.colorbar()
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/librosa/display.py:657: UserWarning: Trying to display complex-valued input. Showing magnitude instead.
  warnings.warn('Trying to display complex-valued input. '
Out[26]:
<matplotlib.colorbar.Colorbar at 0x1159a3128>
In [27]:
# STFTの結果Dからクロマベクトルを算出
chroma_stft = librosa.feature.chroma_stft(S=D, sr=sr)

fig = plt.figure(1, figsize=(12,4))
ax = fig.add_subplot(2,1,1)
librosa.display.specshow(chroma_stft, y_axis='chroma')
plt.title('chroma_stft'); plt.colorbar()
ax= fig.add_subplot(2,1,2)
librosa.display.specshow(chroma, y_axis='chroma', x_axis='time')
plt.title('chroma_cqt'); plt.colorbar()
plt.tight_layout()
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/librosa/display.py:657: UserWarning: Trying to display complex-valued input. Showing magnitude instead.
  warnings.warn('Trying to display complex-valued input. '

その他特徴量抽出

  • メルスペクトル、MFCC、tonnetzなども算出可能
  • $\Delta$特徴量の計算も可能
In [28]:
M = librosa.feature.melspectrogram(y=y, sr=sr)
mfcc = librosa.feature.mfcc(y=y, sr=sr)
tonnetz = librosa.feature.tonnetz(y=y, sr=sr)

librosa.effectsモジュール

  • 音源分離
    • librosa.hpss
      • 調波成分と打音成分に分解
  • 時間伸縮、ピッチシフト
    • librosa.time_stretch
    • librosa.pitch_shift
In [29]:
y_harmonic, y_percussive = librosa.effects.hpss(y)
In [30]:
IPython.display.Audio(y, rate=sr)
Out[30]:
In [31]:
IPython.display.Audio(y_harmonic, rate=sr)
Out[31]:
In [32]:
IPython.display.Audio(y_percussive, rate=sr)
Out[32]:
In [33]:
# それぞれに対するCQTを求める
C      = librosa.cqt(y=y, sr=sr)
C_harm = librosa.cqt(y=y_harmonic, sr=sr)
C_perc = librosa.cqt(y=y_percussive, sr=sr)
In [34]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power= librosa.amplitude_to_db(np.abs(C), ref=np.max)
log_cqt_harm = librosa.amplitude_to_db(np.abs(C_harm), ref=np.max(np.abs(C)))
log_cqt_perc = librosa.amplitude_to_db(np.abs(C_perc), ref=np.max(np.abs(C)))

ax = fig.add_subplot(3,1,1); librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
ax = fig.add_subplot(3,1,2); librosa.display.specshow(log_cqt_harm , x_axis="time", y_axis="cqt_hz")
ax = fig.add_subplot(3,1,3); librosa.display.specshow(log_cqt_perc , x_axis="time", y_axis="cqt_hz")
Out[34]:
<matplotlib.axes._subplots.AxesSubplot at 0x1139ae2e8>
In [35]:
# 2倍の速度に時間伸縮
y_fast = librosa.effects.time_stretch(y, 2.0)
IPython.display.Audio(y_fast, rate=sr)
Out[35]:
In [36]:
# 3度上にピッチシフト
y_third = librosa.effects.pitch_shift(y, sr, n_steps=4)
IPython.display.Audio(y_third, rate=sr)
Out[36]:
In [37]:
# 減五度にピッチシフト
y_tritone = librosa.effects.pitch_shift(y, sr, n_steps=-6)
IPython.display.Audio(y_tritone, rate=sr)
Out[37]:

onset検出とビートトラッキング

  • librosa.onsetモジュール、librosa.beatモジュールを使用
    • 出音の強度をonset.onset_strengthで算出
    • 強度を元にonset箇所をonset.onset_detectで検出
    • 強度を元にbeat.beat_trackでテンポとビートを推定
In [38]:
onset_envelope = librosa.onset.onset_strength(y, sr)
onsets = librosa.onset.onset_detect(onset_envelope=onset_envelope)
In [39]:
fig = plt.figure(1, figsize=(12,8))
ax = fig.add_subplot(2, 1, 1)
ax.plot(onset_envelope, label="onset strength")
ax.vlines(onsets, 0, np.max(onset_envelope), color="red", alpha=0.25, label="onsets")
ax.set_xlim([0, len(onset_envelope)]); ax.legend(loc=1)
ax = fig.add_subplot(2, 1, 2)
librosa.display.waveplot(y, sr)
Out[39]:
<matplotlib.collections.PolyCollection at 0x113aa6c88>
In [40]:
tempo, beats = librosa.beat.beat_track(onset_envelope=onset_envelope)
In [41]:
print(tempo)
129.19921875
In [42]:
fig = plt.figure(1, figsize=(12,8))
ax = fig.add_subplot(1, 1, 1)
ax.plot(onset_envelope, label="onset strength")
ax.vlines(onsets, 0, 0.8 * np.max(onset_envelope), color="red", alpha=0.5, linestyle="dotted", label="onsets")
ax.vlines(beats, 0, np.max(onset_envelope), color="black", alpha=0.25, label="beats")
ax.set_xlim([0, len(onset_envelope)]); ax.legend(loc=1)
Out[42]:
<matplotlib.legend.Legend at 0x114fc1128>

クリック音を生成し同時に再生

  • MIRの評価に用いられるmir_evalライブラリを使用
In [43]:
!pip install mir_eval
Requirement already satisfied: mir_eval in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (0.4)
Requirement already satisfied: six in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.11.0)
Requirement already satisfied: numpy>=1.7.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.14.2)
Requirement already satisfied: scipy>=0.9.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.1.0)
Requirement already satisfied: future in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (0.16.0)
In [44]:
# http://colinraffel.com/publications/ismir2014mir_eval.pdf
import mir_eval
In [45]:
beat_times = librosa.frames_to_time(beats)
y_click = mir_eval.sonify.clicks(beat_times, sr, length=len(y))
In [46]:
IPython.display.Audio(y_click, rate=sr)
Out[46]:
In [47]:
IPython.display.Audio(y + y_click, rate=sr)
Out[47]:

特徴量のビート同期

  • librosa.util.syncを用いるとビートとその他特徴量を同期させることができる
    • 各ビートの間の成分を平均化したり中央値を取ったりできる
In [48]:
tempo, beat_frames = librosa.beat.beat_track(y=y_percussive, sr=sr)
# 13次元のMFCCおよびΔを抽出
mfcc = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)
mfcc_delta = librosa.feature.delta(mfcc)
In [49]:
print(beat_frames.shape)
print(mfcc.shape)
print(mfcc_delta.shape)
(98,)
(13, 2647)
(13, 2647)
In [50]:
# ビート間の平均値を用いてMFCCとΔMFCCをビートに同期
beat_mfcc_delta = librosa.util.sync(np.vstack([mfcc, mfcc_delta]),
                                    beat_frames)
In [51]:
beat_mfcc_delta.shape
Out[51]:
(26, 99)
In [52]:
# 調波成分からクロマベクトルを推定
chromagram = librosa.feature.chroma_cqt(y=y_harmonic, sr=sr)

# ビート間の中央値を用いてクロマベクトルをビートに同期
beat_chroma = librosa.util.sync(chromagram, beat_frames, aggregate=np.median)

# ビート同期特徴量(クロマベクトル、MFCC、ΔMFCC)を結合
beat_features = np.vstack([beat_chroma, beat_mfcc_delta])
In [53]:
# MFCC 13dim, ΔMFCC 13dim, クロマベクトル(12音階; 12dim)
beat_features.shape
Out[53]:
(38, 99)
In [54]:
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(211)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
log_beat_features = librosa.amplitude_to_db(np.abs(beat_features), ref=np.max)
ax = fig.add_subplot(212)
ax.pcolormesh(log_beat_features)
Out[54]:
<matplotlib.collections.QuadMesh at 0x113a0ac50>

まとめ

  • librosaを用いた音楽情報検索(MIR)について説明
    • 特徴量抽出、可視化、音源分離、ビート推定
  • 音楽データに対する探索的分析の手段の一つ

以下補足およびおまけ

librosa.decomposeによる音源分離

  • NMFによるモノラル音源分離を行う
    • スペクトログラムを行列と見てscikit-learnのNMFに投げる
    • アクティベーション行列と基底行列に分解
    • 教師なし音源分離が可能

NMF

  • 非負値行列$\textbf{X}$を$\textbf{X} \simeq \textbf{WH} $と分解する手法
    • Hを基底行列、Wをアクティベーションと呼ぶ
  • scikit-learnのNMFでは下記目的関数を最小化

$$0.5 * ||\textbf{X} - \textbf{WH}||_{Fro}^2

  • \alpha l1_{ratio} ||vec(W)||_1
  • \alpha l1_{ratio} ||vec(H)||_1
  • 0.5 \alpha (1 - l1{ratio}) ||W||{Fro}^2
  • 0.5 \alpha (1 - l1{ratio}) ||H||{Fro}^2 $$

  • 通常最適化には補助関数法が使われる

In [55]:
D = librosa.stft(y)
S, phase = librosa.magphase(D)
comps, acts = librosa.decompose.decompose(S, n_components=8, sort=True)
S_approx = comps @ acts
In [56]:
print(comps.shape)
print(acts.shape)
print(S.shape, S_approx.shape)
(1025, 8)
(8, 2647)
(1025, 2647) (1025, 2647)
In [57]:
def plot_decomposed(S, comps, acts, S_approx):
    fig = plt.figure(figsize=(10,8))
    ax = fig.add_subplot(311)
    librosa.display.specshow(librosa.amplitude_to_db(S, ref=np.max), y_axis='log', x_axis='time')
    ax.set_title('Input spectrogram'); plt.colorbar(format='%+2.0f dB')

    ax = fig.add_subplot(323)
    librosa.display.specshow(librosa.amplitude_to_db(comps, ref=np.max), y_axis='log')
    plt.colorbar(format='%+2.0f dB'); ax.set_title('Components')

    ax = fig.add_subplot(324)
    librosa.display.specshow(acts, x_axis='time')
    ax.set_ylabel('Components'); ax.set_title('Activations')
    plt.colorbar()

    ax = fig.add_subplot(313)
    librosa.display.specshow(librosa.amplitude_to_db(S_approx, ref=np.max), y_axis='log', x_axis='time')
    plt.colorbar(format='%+2.0f dB')
    ax.set_title('Reconstructed spectrogram')

    fig.tight_layout()
In [58]:
plot_decomposed(S, comps, acts, S_approx)
In [59]:
S_partly = comps[:, :3] @ acts[:3, :]
In [60]:
plot_decomposed(S, comps[:, :3], acts[:3, :], S_partly)
In [61]:
y_partly = librosa.istft(S_partly * phase)
IPython.display.Audio(y_partly, rate=sr)
Out[61]:

decomposeの用途

  • 視聴用途の音源分離には(このままでは)向かない
    • 振幅スペクトルのみを分解するが、位相情報は元信号の元を使うため歪が大きい
    • 基底数を増加した上で適切なクラスタリングを施して使用したり、事前学習した楽器別の基底などを用いて分離を行うなどの工夫が必要
  • 音響イベント検出のための特徴量に用いることも可能